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Abstract 

We introduce a new probabilistic method for solving a class of impulse control 
problems based on their representations as Backward Stochastic Differential Equations 
(BSDEs for short) with constrained jumps. As an example, our method is used for pric- 
ing Swing options. We deal with the jump constraint by a penalization procedure and 
apply a discrete-time backward scheme to the resulting penalized BSDE with jumps. 
We study the convergence of this numerical method, with respect to the main approx- 
imation parameters: the jump intensity A, the penalization parameter p > and the 
time step. In particular, we obtain a convergence rate of the error due to penalization 
of order (Ap)"~2,Va G (O, ^). Combining this approach with Monte Carlo techniques, 
we then work out the valuation problem of (normalized) Swing options in the Black 
and Scholes framework. We present numerical tests and compare our results with a 
classical iteration method. 
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1 Introduction 



In this report, we introduce a new probabilistic method for solving impulse control 
problems based on their representations as Backward Stochastic Differential Equations 
(BSDEs for short) with constrained jumps. As an example, our method is used for 
pricing Swing options in the Black and Scholes framework. 

BSDEs provide alternative characterizations of the solution to multiple-obstacle, 
optimal switching (see among others [HI HI [20l [26l [13]) and more generally impulse 
control problems: Kharroubi et al. [21] recently introduced a family of BSDEs with 
constrained jumps providing a representation of the solution to such problems. A 
challenging question is that of the numerical approximation of this kind of BSDEs 
with constrained jumps. 

A discrete-time backward scheme for solving BSDEs with jumps (without con- 
straint) has been introduced by Bouchard and Elie j3j. In our case, the main difficulty 
comes from the constraint, which concerns the jump component of the solution. These 
BSDEs do not a priori involve any Skorohod type minimality condition. In consequence, 
classical approaches by projected schemes (discretely reflected backward schemes) used 
for example by [2] and [TU] cannot be used. 

We consider a penalization procedure to deal with the constraint on jumps and 
provide a convergence rate of the penalized solution to the exact solution. This allows 
us to establish a convergence rate of the error between the solution of the considered 
impulse control problem and the numerical approximation given by the discrete-time 
solution to the penalized BSDE with jumps, as the penalization coefficient and the 
number of time steps go to infinity. 

The rest of the report is structured as follows: in Section [21 we set the considered 
impulse control problem in the mathematical framework of BSDEs with constrained 
jumps. We present in Section [3] our penalization approach and provide a global con- 
vergence rate of our approximation. In Section |31 our method is used for pricing 
multi-exercise options, so-called (normalized) Swing options. This multiple optimal 
stopping time problem leads to a particularly degenerate three-dimensional impulse 
control problem. We combine our BSDE-based approach with Monte Carlo techniques 
and deal with Swing options with a small maximal number of exercises rights, due to 
large computational times. We compare our pricing results with those obtained by a 
classical iteration-based approach proposed for example by [9] . 

2 BSDE Representation for Impulse Control Prob- 
lems 

Let T be a given time horizon. We work in a complete probability space (51, J-", P), on 
which is defined a c?-dimensional Brownian motion W and a Poisson process N with 
intensity A > 0. We denote by F = {J^t)t>Oi the augmentation of the natural filtration 
generated by W and iV, by = (J^('^)t>o the one generated by W , and by V, the 
cr- algebra of predictable sub-sets of 17 x [0,r]. 
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Notation 



Throughout this report, the eucUdean norm defined on M.'^ or on M will be indis- 
criminately denoted by I • |. The matrix transposition is denoted by _L. In addition, 
unless specified otherwise, C will denote a strictly positive constant depending only on 
Lipschitz constants of the coefficients of the problem, see assumptions (H) and (H') 
below, and constants T, |5(0)|, |cr(0)|, |7(0)|, |/(0)|, |«;(0)| and \g{0)\. 

Besides, we shall use the standard notations: 

• S^, the set of real- valued cadlag adapted processes Y — {Yt)o<t<T such that 



IFII 



E 



sup \Yt\ 

0<t<T 



< 00 



• A"^, the sub-set of <S^ such that 



A"^ := {K e 5^ : {Kt)Q<L<T nondecreasing , Kq = 0} . 
• L|([0,T]), the set of real-valued adapted processes {4>t)o<t<T such that 



E 



\(l>t\^dt 



< 00 . 



L'^{W), the set of real-valued 'P-measurable processes Z = {Zt)s<t<r such that 









/ \Ztfdt 




Jo 



< 00 . 



L'^{N), the set of real-valued 'P-measurable processes V = {Vt)s<t<r such that 



Jo 



< 00 . 



• V denotes the set of P-measurable essentially bounded processes, valued in (0, oo) 
and VP = {i/P G V : < p a.s.}. 

2.1 A Class of Impulse Control Problems 

We consider the class of impulse control problems whose value function is defined by: 



v{t,x) = sup E 

«=('r)=));>l6^t,T] 



0+ r f{Xl'=^nds+ Yl 
Jt 1~1 



k>l 
t<Tk<T 



(1) 
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An impulse strategy u — {Tk)k>i is said to be admissible for problem ([Ij (and belongs 
to U(t^T]) if it is a non-decreasing sequence of F^-stopping times valued in {t,T] (we 
set by convention tq = t) such that, if 

n^t,T] ■■=Hk>l-t<Tk<T} 

denotes the (random) number of interventions of the strategy u before time T, then 



E 



"•(t,T] 



2 



< C, (2) 



for some universal constant C > 0. The controlled state variable X*'^'" is a cadlag 
process with dynamics 



r 6(X*^"'")dr + r a(X*^"'")dW, + V 7(^*'-'"), Vs > t. (3) 



Between two successive intervention times and Tk+i, the state variable evolves as a 
diffusion process and the controller makes an integral profit /. At each decided inter- 
vention time Tk, he gives an impulse to the system: the state process jumps with a size 
X" — = 7(Ar"_) and he obtains the intervention gain k. 

We consider standard assumptions on the coefficients of the problem: 

(H) 6 : R'^ 1-^. R'^, cr : E'^ H> K''^'^ and 7 : R'* M- R** are Lipschitz continuous 
and 7 is uniformly bounded. 

/ : R'' R, K : R"* R and .g : R'' R are Lipschitz continous. 

(H') The maps 6, fi, 7, and g belong to C^iW^) and have Lipschitz continuous 
derivatives. 

A straightforward computation using (H), ([2]) and Gronwall's lemma shows that 



V(<, x) e [0, T] X R'*, sup E 
u£U^tT] 



sup \X 

t<s<T 



t'XAi\ 



< 00. (4) 



Finally, we will assume the existence of an optimal strategy u* — {T^)k>i G l^(t.T\ 
to problem ([T|). We refer for example to 4J and ^25] in the infinite horizon case, for 
specific conditions on the coefficients of the problem which ensures such an existence. 

2.2 Link to BSDEs with Constrained Jumps 

Let us consider the BSDE with constrained jumps 

{Yt^ giXr) + f{X,)ds - ZsdWs - V^dNs -f dK^, WO<t<T 
\vt + K{Xt-) <0, VO<i<T 

(5) 

where X is the (uncontrolled) jump diffusion process with dynamics 

dXt = biXt)dt + a{Xt)dWt + -f{Xt-)dNt. (6) 
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Under (H) , this SDE admits an unique solution in and it is shown in |21| that under 
the additional assumption (Hi) given below, ([5]) admits a unique minimal solution 
(Y, Z, V, K)gS^ X h^{W) X L2(iV) X ^2 ^ith K predictable. 

The solution {Y, Z, V, K) is said to be minimal if and only if it has the smallest 
component Y in the (infinite) class of solutions to ([S]). {Yt)t>Q is called the value 
process and jumps with a size Vt — Yt — Yi- . 

(Hi) There exists a solution [Y , Z, K) e S'^ x L'^{W) x to 

Yt=g{XT) + f{X,)ds- ZsdWs+J^ K{X,- )dNs+J^ dK,. 
(H'l) (Hi) holds and Yt — v{t, Xj), VO < < < T for some v with linear growth. 

(H2) There exists a non negative function (p e C^(R'') and a constant p > s.t. 

^'P + / < Pf, 95 ^ T~(-(p > 0, 

<y5 > ff, lini|2;|^oo = 00, 

in which C is the local component of the generator of the process X and H, the 
intervention operator: 

jCv{t, x) = b{x) ■ D^v{t, x) + ^Tr {c7C7^{x)Dlv{t, x)) , 
'Hv{t, x) = v{t, X + 7(2;)) + k{x) . 

Let (y,*'=^,Z*^^,V;*^^,i^*'^)t<^<T be the solution to © when X = (X*^^)t<,<T is 
the solution starting at a; in i to SDE Under assumptions (H), (H'l) and (H2), [3T] 
show that the solution to impulse control problem ([!]) coincides with initial value of 
component y*'^: 

y/- = t.(i,x) (7) 

and is equal to the (unique) solution with linear growth to quasi- variationnal inequality 

'^^{"IK^'^) ~ ^v{t,x) - f{t,x) ; 

v{t,x) -nv{t,x)} = 0, V(t,x) e [0,T) X R'*, (8) 
min{v{T-,x) - g{x);v{T-,x) -nv{T-,x)} = 0, Vx G R'^. 

Let us mention that in general, the terminal condition v{T^,-) — g is irrelevant, 
because of the possible discontinuity of Y in T~ due to constraints: the relaxed terminal 
condition in 1^ expresses the possibility of a jump at time . 

Remark 1. For a better intuition, the following interpretation to solution (Y, Z, V, K) 
holds when assuming u G C^'^([0, T], R''): 

VO<t<r, Yt = v{t,Xt) 

Zt = a{t,Xt-)D.,v{t,Xt-) 

Vt = v{t,Xt- +^{Xt-))-v{t,Xt-) 

= nv{t,Xt-)-v{t,Xt-)-K{Xt-) 
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The constraint in ([5]) means thus that the obstacle condition is satisfied, namely 

v{t,Xt-)-nv{t,Xt-)>o. 



3 Convergence of the Numerical Approximation by 
Penalization 

It does not seem possible to use the minimality condition of the solution to BSDE with 
constrained jumps ([5]) directly in a numerical scheme. We thus propose an approach 
by penalization of the jump constraint. The penalized constraint is introduced in the 
BSDE driver: when the constraint is fulfilled, this penalization term disappears, and 
otherwise penalizes the driver with an exploding factor p. 

In Theorem [U we provide an explicit rate of convergence of our approximation 
with respect to the parameters introduced: namely, the jump intensity A, the penaliza- 
tion coefficient p and the time step. Such an error estimate is essential for numerical 
purposes (understanding of the numerical impact of those parameters) and allows to 
adjust in practice the fineness of the time grid in relation to (A,p). 

3.1 Approximation by Penalization 

Given a parameter value p > Q, the penalized BSDE is: 

Yr = g{XT) + j"^ [f{Xs) + p {VP + k{X,- ))+ a] ds (9) 
ZfdWs-J^ VPdN,, VO<i<T 

which admits an unique solution {YP,ZP,VP) G 5^ x L^{W) x L^(iV) from the clas- 
sical theory of BSDEs with jumps. In addition, the sequence of penalized solutions 
(YP, ZP, VP)p tends in L2([0, T]) x L^{W) x L^{N) to the minimal solution (Y, Z, V) to 
([5]) as p goes to infinity, see [21]. Besides, the convergence of {YP)p to Y is monotone 
and increasing. 

Let {YP^t,x^2p,ux^yp,t,x^^ ^^le solution to ^ when X = (X*^^)t<^<T. We con- 
sider the following error introduced by this penalization procedure: 

£P := sup (10) 

0<t<T 



For any i < 77 < T, let us introduce 



w^(i, x) — sup E 

"=(Tfc)fc>ieW(t_T-,] 



k>\ 
t<Tk<T 



(11) 



which corresponds to initial problem ([T]) restricted to the sub-set of strategies taking 
values in {t,T — 77]. Wc shall denote by u^* = {T^*)k>i an 772-optimal strategy to 
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problem pT|) (the existence of an optimal strategy is not ensured) and by n^* be the 
number of impulses in strategy u^* that is: 

n"* := tt {fc > 1 : i < rf < T - 77} . 

We will use the following additional assumptions: 
(H") There exists some n e N* such that 

Vi>n, PK* > j) <;(j) 

for some map I such that l{j) < e'^^^ for some some constant C > 0. 
(H*) There exists a map h such that h{e) = O^^oie^) and 

Ve > 0, P (min \r^,l, - | < < h{e). 

Remark 2 (Assumptions (H") and (H*)). Both assumptions (H") and (H*) are di- 
rectly satisfied for the problem of Swing options valuation since the number of exercises 
right is almost surely bounded by some rimax and there is some fixed time delay 6 > 
between two consecutive interventions. 

More generally, (H") is intuitively satisfied as soon as the controlled state variable 
is constrained almost surely and admits jumps of constant sign, see Example [1] (H*) 
is verified for sufficiently smooth problems, see for example the case of optimal forest 
management studied in [Tl. 

Example 1. Let us assume that the state variable defined in ([3]) is such that 

• 6 is uniformly bounded and ct > constant, 

• for some constant c > 0, 

sup 7(2:) < — c, 

and that the optimal strategy u* implies > a.s. Then a straightforward com- 
putation shows that the (random) number rt^g j,j of optimal impulses before time T 
satisfies, for any a > 0, 

P > = O (e"°") as n ^ +00. 



Proposition 1. Let assumptions (H), (H"), (H*), (H'l) and (H2) be satisfied. Then 
the penalization error in (jTU]) admits the following bound as p goes to infinity: 

for some constants C > and C > 1, which do not depend either on X, p, n or a. 
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Proof. We provide the main arguments of the proof and refer the reader to [1_ for more 
details. The main idea comes from the following exphcit fmictional representation 
available for FP-*-^, see [2T|: 



ess supl 



(12) 



where E"'' denotes the expectation under the probability measure P'^'' equivalent to 
on (51,7^7^) with Radon-Nikodym density 



dF 



-/o^K-l)Ads„/o^lnK)dAr, 



The specificity of such a change of measure is that it impacts only the jump parts of the 
processes: under P" , the Brownian motion W remains unchanged whereas N has a 
(stochastic) intensity (Ai'f )s>o- We shall denote by N^, the doubly stochastic Poisson 
process (Cox process) with intensity (Az^f )s>o under P, by (r^)fe>i the sequence of its 
jump dates and by the solution to © driven by N^. 

In view of (jlip and p2p . we introduce a convenient measure change, which intu- 
itively forces the penalized solution to jump as soon as possible after that an optimal 
impulse occurs: 

1 else. 



Vs>0, lyP 



(13) 



Notice that i/^ is a "P-measurable process bounded by p a.s. By definition of the 
counting process iV^, 



Vs > 0, 



In other words, the increment — r^* has an exponential distribution with parameter 
(Ap), conditionally to the fact that X^ has jumped one time less than X"'' . This 
allows us to compute an estimate of the distance 

sup |w^(t,a;) - , 

0<t<T 



which holds for by sending rj to together with a continuity argument of the value 
function in its maturity variable. □ □ 



3.2 Convergence Rate of the Numerical Scheme 

Given a regular time grid tt ~ {to = 0, ii, . . . , t^v = T}, we assume that the solution X 
to ([6]) can be simulated on tt either perfectly or by using an Euler scheme and denote its 
discrete-time version by X'^. Along the lines of [3], we consider the following backward 
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discrete-time scheme for numerically solving the penalized BSDE with jumps ([S]), 



yP:^ _ 1 



- Ai„+i 



AAt„ + i 

1 E 



E[y,^';^A7v,„,jj-,, 



(14) 



At 



n+1 



where Atn+i = tn+i — tm AWt^^-^ is the Brownian increment on and ANt^_^_^ 

the compensated version of the Poisson increment ANt^^-^ on [i„,t„+i). 

We consider the classical discretization error between the continuous-time solution 
(FP, ZP, VP) in ^ and its discrete-time approximation (FP-'^, Z^^'^, yP''^) in HI]), that 
is 



£''(yP) := max sup E\YP~YP-''\ 



r{ZP) := J2 



^N-l „t„ + i 



^n=0 



£^{VP) := \^J2J IE I - F/;" I ' Adi j . 

Because of the lack of first order regularity of the driver of the penalized BSDE 

fP{x,v):^ f{x) + {p{v + K{x))+ -v)X, V(x, v) G M"* X M, (15) 

classical regularization arguments for the FBSDE coefficients and Malliavin differenti- 
ation representations allow us to provide an explicit convergence rate of order |7r| = for 
errors £'^{Yp) and £'^{Vp), but only of order |7r|3 for error £'^{Zp), see Proposition [2] 

The impact of the penalization coefficient p on the convergence of backward discrete- 
time schemes is well known in practice, even if there does not exist, to our best knowl- 
edge, any explicit computation in the literature (see for example the numerical experi- 
ments of [2^ for the resolution by penalization of a BSDE with one reflecting barrier). 
Basically, as p increases at fixed discrete-time step, the quantity fP{-)Atn+i explodes, 
leading to a numerical explosion of the approximate values YP''^, see (|14p. 

We show rigorously in Proposition [5] that the discretization error grows exponen- 
tially with (Ap^). This is due to the linear dependence in A and p of the BSDE driver 
fP, see (|15l) . and estimate computations based on the use of Gronwall's lemma. 



Proposition 2. At 



(H). Then, as soon as 



1 



(16) 
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we get the following bounds as p goes to infinity: 



S^'iVP) < c((l + A)2A3p2(7Ap^|^|f 
Under (H'), there exists a version of such that, 

for some constants C > and C > \, which do not depend either on \, p or |7r|. 

Proof. This follows from the same arguments as |15j : computations using ltd and 
Gronwall's lemma and regularization and Malliavin differentiation arguments applied 
to the penalized BSDE with jumps ([9]). We refer the reader to |lj for a detailed 
proof. □ □ 

Propositions [T] and [2] enable us to establish a global convergence rate of the error 
introduced by our approximation by penalization. 

Theorem 1. Let the assumptions of Proposition]^ be satisfied. Then 

EP + £^YP)<cl^j^^^ + {l + \f\pC>^p"\^\^, VaG(^O,0 (17) 

for some constants C > and C > 1, which do not depend either on \, p, \t:\, n or a. 
Thus, for a sufficiently small time step |7r| with respect to A and p, the global error is 
such that 



Remark 3 (Global convergence rate). At fixed time step |7r|, the convergence rate 
strongly deteriorates as A or p increases, see (|17p . The numerical method is more 
sensitive to p than to A according to and ^T7\ . constitutes a necessary 

condition for the convergence of the backward discrete-time scheme. In practice, the 
penalization parameter will need to be chosen relatively small and the time step |7r| 
very small to avoid multiple jump times on each time step (otherwise, this introduces 
a bias). 



4 Application to Swing Options Valuation 

In this section, our method is applied for the valuation of Swing options in the Black and 
Scholes framework. This constitutes a multiple optimal stopping time problem which 
can be reformulated as a particularly degenerate three-dimensional impulse control 
problem. We have been able to achieve convergence for a small maximal number of 
exercises rights (jT-max < 2) due to the slow computational speed of our method. 
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4.1 Swing Options Valuation as an Impulse Control Problem 

We consider a (normalized) Swing option: the holder of the option is given a maximal 
number of exercise rights, say n^ax ^ 1; a-nd has the opportunity to sell whenever he 
wants over a time period [0, T] an underlying asset against a fixed strike price. 

For some fixed strike price K, we shall denote by <j){s) = [K — s)+ the reward 
function corresponding to the profit made at each exercise date and by 5* the underlying 
asset spot price. We concentrate here on the risk-neutral Black and Scholes framework 
in which r > is a constant interest rate and the spot price process is defined by 



St = s exp <; (r - ^CT^)* + aWt 



yt > 0, 



(18) 



where a > denotes the volatility coefficient and s an initial price value. 



A delay 6 > between two consecutive exercise dates is introduced. Indeed, without 
any delay, the optimal strategy would consist in rimax simultaneous exercise at an 
imique optimal date (so that this option is equivalent to rimax identical American 
options). The value of such an option at time can be written as the solution to the 
following multiple optimal stopping time problem 



sup 
M=(Tfc)fc>iew 



E 



(0,T] 



fc>i 



(19) 



in which a strategy u = {Tk)k>i is said to be admissible and belongs to U^^ j,j if and 
only if it is a vector of F^-stopping times valued in (0, T] with maximal lenght rimax 
which satisfied the constraint on delay, that is 



Vfc > 1, Tfe4 



Tk > S. 



As a multiple optimal stopping time problem, this problem can be formulated as 
a particular (and strongly degenerate) impulse control problem. An impulse control 
corresponds to a sequence of exercise dates and the intervention gain is written as 
the payoff function <p multiplied by an indicator function which allows to satisfy both 
the constraint on the number of exercise rights and the constraint on delay between 
exercise dates. Namely: 



v{s) 



Slip 

"=(Tfc)fc>i 



J2 e-'-^''</'(5.Jl 



k>l 
■Tk<T 



{(e^:^>5)n(Q^_<n„..)} 



(20) 



where u = {Tk)k>i is the sequence of F'^-stopping times valued in (0,r] and two ad- 
ditional state variables which both are controlled and discontinuous (cadlag) processes 
are introduced: 

• counts the number of exercise rights used before considered time 

(3o=0, Qt =Hk>l,rk<t}, Vt>0. 
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• 0" := — inf {t — Tk,Tk < t} corresponds to the delay between t and last 
exercise date 



where by convention 8q = 9q = 0. 
Obviously, the problem is equivalent to the impulse control problem (1^0)) : 

VseM, w("""'-)(s) =w(s). 

Penalized BSDE Associated to Swing Option Pricing Problem 

Let us introduce the uncontrolled variable {Q, Q) defined by 

iQt = Nt, yt>o, 



(21) 



where we recall that iV is a Poisson process with intensity A > 0. For a penalization 
coefficient p > 0, the penalized BSDE with jumps associated to problem (PO)) is 



Z^dWu 



VO<t<T. 



(22) 



in which the intervention gain is defined by 

K(s,q,0) :==0(s)l{(e>5)n(g<„_.-i)}, V(s, ^?) £ M x N X R+. (23) 

4.2 Numerical Valuation Algorithm 

The three-dimensional process {S, Q, Q) can be exactly computed on the time grid 
TT. We shall denote by {S"^ its version on tt. In particular, the pure jump 
processes Q and O can be computed without approximation error, recall (j21l) . from a 
simulated trajectory of a Poisson process with intensity A, see for example [11 . 

As the driver of the penalized BSDE with jumps (p2)) does not depend on Z^, it is 
sufficient to compute {YP-^ , V^''^) on tt, which can be made by the backward recursive 
scheme: 

yt„ G TT, t„ < T : 



1 



t„ l+rAt„+i 



{e [Yr::,\Tt„' 



(24) 



AAt 



n+l 
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and in our setting: 

E[-\Tt„]^E[-\iSl,Ql,ei)] 
Monte Carlo-based Resolution 



(25) 



We compute estimators of the conditional expectations by a classical least squares 
Monte Carlo technique. We use least squares regressions on adaptative local basis 
functions, see [5]. Since such an approach is only relevant for real- valued variables (the 
regression basis functions have compact support), it cannot handle the integer- valued 
variable Q'^. 

We thus deal with this variable explicitly. Namely, we simulate M > 1 i.i.d. paths 

of (S"^,g'^,e'^) 

(5,7r,(m)^ g7r,(m)^ Q7r,(m))^ Vm < M, 

and at each time step i„ < T of the backward recursion, the set of M Monte Carlo 
samples is separated in nmax + 1 sub-sets corresponding to the samples on which Q^^ = 
0, 1, . . . , rimax-i and Q^^ > nmax- Then, we just need to estimate rimax conditional 
expectation operators, namely 

E [• I {SI , Ql - g, erj] , Vg e {0, 1, . . . , n,nax - 1} 

by using the corresponding Monte Carlo samples, since {Yf^^ ,Vf^^) — (0,0) when 
Qtn — ^max, see Remark 21 

Remark 4. Let us show that > rimax implies Yf''^ = Vf'^ = 0. We have 



E 

= E 



y/^;^ ANt„^, I {SI , Ql > n„,ax, ei 



E 



Yr:^,\Ql > n„,ax ANt^^,\{Sl,Ql > n^ax,e7J 



0, 



by definition of k in (|23p and ([24]) . In the same way 



E 



which allows to conlude. 

The Monte Carlo-based algorithm is then the following 
I. Initialization; 

r/„-('"' = <^(5-('"))i 



0}' 



Vto < M. 



II. Computation backward in time of (V^P''^'('")j yP'^'(™)) on each sample m < M. 
For n = TV - 1, ... ,0, set 



Ml TO=1, 



,M:Ql^"'^ 



,M : Q 



max J* ■ 



1, 



Then: 
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1. For any m G Ml 



2. Set q := n,nax - 1- 



V 



,7r,(m) -t^p,7r,(m) 

— ^ + 



= 0. 



3. li q> 0, for any m G Ml , the conditional expectations estimators 



-V,<J,(m) 



Y..q.im) 



E 



E 



are approximated by least squares regression of (vf^^^^^^'' ANj:-'^^_^ 



meMj 



and 



lyP^^rX™) j respectively on 



with b := bl basis functions {f/'ij • ■ • , V'h}- Then, 



AAt 



1 ^"I'.'JX™) 



p,7r,{m) 



p,7r,(m) 



AAt, 



n+l 



4. g := g — 1 and go to 3. 

III. At time to {M^to = {1, ■ • ■ , ^^1 and Ml^ 
estimator is given by Y^^^ such that 



, > 1) the Swing option price 



^to — AAti Af Z^m=l y^ti ^^^ii 

Let us highlight some features of the above-presented Monte Carlo procedure. At 
each backward induction date i„ < T, we have to estimate in worst cases 2 X ^T-niax 
conditional expectations, performed on each subset M-1^,q < n^nax — 1- When n„iax 
increases, much more Monte Carlo samples are needed as each least squares regression 
requires a sufficient number of samples. 

In addition, the number of local basis functions bf^ has to be adapted to the number 
of Monte Carlo samples used for the least squares regression, namely card(A4j^). We 
thus introduce a dynamic choice for it is fixed proportionally to card(A^^^) for 
any q < n,„ax - 1 and t„ < T. 

Remark 5 (Statistical error of our method). A control of the statistical error intro- 
duced by the least squares Monte Carlo approach is provided in Gobet et al. [TH] (see 
[33] for further details). By extension, this applies to BSDE with jumps (see [15]) to 
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control the error on the jump component y^'^ and thus ensures that the least squares 
Monte Carlo error tends to as the number of samples M and the number of basis 
functions h tends to +00. 

A Benchmark Method Based on Iteration 



The classical method to value such a Swing option, recall formulation HH), is based 
on an iteration over the number of exercise rights, see for example [9]. The dynamic 
programming principle provides a direct link between the solution v^^^ to the same 
problem as pop but with at most j < rimax exercise rights and the solution w^^"^) with 
at most (j — 1) exercise rights. 

The value of the Swing option with exercise right is obviously zero w^"-* = 0. 
Then, we compute the sequence of values of Swing options with j exercise rights u*^^^ , 
j = 1, . . . , ninax according the backward recursion scheme: 



Vt„ G 7r,T-(5< t„ < T : 

v^^Htn,s) =max|(/i(s) ; e-'''^*"+iE(*-'*) 

\ftn G Tr,tn < T ~ S : 

?;W(i„,s) =max{0(s)+e-^^E(*"'^) [v^^-^^t,, + S, Sl^g)] 



rAt„ + i|p(t„,s) 



(^n+li^T^+i) } 



where E^*"'*' [■] := E [■|5'J^ = s] • We use the same least squares Monte Carlo regression- 
based method for approximating the conditional expectations operators as above. 



4.3 Pricing Results 

We consider put options with maturity T = 1 year and a strike price K = 100. The 
Black and Scholes parameters, see are r = 0.05, a = 0.3 and s = 100. 



4.3.1 Special Case of American Options: nmax = 1 

In the single-exercise case, the additional variable Q disappears (there is no delay con- 
straint). This helps simplify the Monte Carlo procedure described in Paragraph 14.21 
In particular, the algorithm implies only two sequences of samples subsets, whether 
one jump of N occurs before T on the considered path or not: that is ((7W^)t„)ri<Af-i 
and ((A^")t„)„<Ar_i. 

In our numerical experiments, we find out that increasing too much A makes the 
variance of the Monte Carlo procedure explode. It would be necessary to increase 
the number of Monte Carlo samples, which leads to prohibitive computational times 
(each pricing result presented below was obtained after a computation between 6 and 
8 hours). For the same reason (exploding behavior of the penalized BSDE driver), we 
restrict our numerical experiments to penalization parameters < 5. 
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A 


\ 20 


40 


80 


160 


320 


3 


9.89 


9.92 


9.95 


9.94 


9.83 


4 


9.92 


9.96 


9.99 


9.97 


9.83 


5 


9.95 


9.99 


10.02 


9.98 


9.76 



Table 1: Approximate prices of an American option with p = 5 



The benchmark price for the American put option is 9.88 (by a binomial approach 
or classical Monte Carlo). We report in Table [T] the price given by our method when 
varying A and the number of time steps N for a penalization parameter equal to 5. We 
used 20 million of Monte Carlo paths. 

In all the experiments that we performed in this simple case, we numerically ob- 
served that the limiting prices of our method (with respect to N) are below the bench- 
mark value: this is due to penalization. 

4.3.2 Swing Options with nmax = 2 
We consider time delays ^ = Jq, j^, Jq- 

The benchmark prices for the Swing put option with 2 exercise rights are 19.27, 
18.77 and 18.21 respectively (computed with the method described in Paragraph l4.21 
N = 200 time steps and Af = 5 million of Monte Carlo paths). We report in Figures 1, 2 
and 3 the corresponding approximate prices when varying A and N for a penalization 
parameter equal to 5 and 10 (we used 40 million of Monte Carlo paths and N = 
20,40,80,160,320,640). 







Bench mark 
Approx. with lambda - 3 
Approx. with lambda - 4 
Approx. with lambda = 5 













































































Benchmark 
Approx. with lambda = 3 
Approx. with lambda - 4 
Approx. with lambda = 5 







































































I 00 200 300 400 500 600 700 

N. number ol" time steps 
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N, number of time steps 



Figm'e 1: Approximate prices of a Swing option with 2 exercise rights 
and S = jQ, with p = 5 (left) and p = 10 (right) 
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Figure 2: Approximate prices of a Swing option with 2 exercise rights 
and S = with p = 5 (left) and p = 10 (right) 
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Benchmark 
Approx. with lambda = 3 
Approx. with lambda = 4 
Approx. with lambda - 5 
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Figure 3: Approximate prices of a Swing option with 2 exercise rights 
and S = with p = 5 (left) and p = 10 (right) 



For each considered value of A, we retrieve a convergence in the number of time 
steps N of our method. As p = 5, approximate prices are converged from TV = 160, 
so that we restrict ourselves to < 320 time steps as p = 10. The limiting values are 
still below the benchmark but accurate option prices (relative error less than 1%) are 
obtained with a penalization coefficient p equal to 10 and N ~ 160. See also Tabled 
in which the (signed) relative error to the benchmark is given in brackets. Besides, we 
observe a monotone convergence in A of our approximate method. 

We should point out that fine-tuning the parameters of the algorithm is difficult. 
As already mentioned, since the number of Monte Carlo paths is different in each set of 
sample paths Mt,^ 1 9 — 0, 1, 2, the number of basis functions used for the least squares 
regressions has to be dynamically adapted. And when increasing much more the jump 
intensity A, more Monte Carlo samples would be necessary. 

For such a Swing option, the running time is much longer because the conditional 
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P 




4 


5 


5 = 


1 

10 


5 

10 


18.80 (-2.44%) 
19.16 (-0.57%) 


18.95 (-1.66%) 
19.24 (-0.16%) 


19.00 (-1.40%) 
19.27 (0.00%) 


5 = 


2 

10 


5 

10 


18.23 (-2.88%) 
18.59 (-0.96%) 


18.37 (-2.13%) 
18.65 (-0.64%) 


18.43 (-1.81%) 
18.69 (-0.43%) 


6 = 


3 
10 


5 

10 


17.62 (-2.97%) 
17.99 (-0.94%) 


17.76 (-2.20%) 
18.04 (-0.66%) 


17.81 (-1.93%) 
18.05 (-0.61%) 



Table 2: Prices of a Swing option with 2 exercise rights 
(hmiting values with N = 160) 



expectations are computed by regression with respect to the bidimensional state vari- 
able (S"^,9'^). The computation of one option price takes at least 15 hours in above 
cases (when N > 80). In comparison, the benchmark method takes less than 5 min- 
utes. Besides, the complexity of our method increases with nmax, leading to untractable 
computational times for bigger values of nmax, see Remark [6l 

On this particular case of Swing options valuation, it seems that our method is less 
competitive than the classical approach. This is without any doubt due to the strong 
degeneracy of such a problem in our impulse control context: the valuation problem is 
3-dimensional and involves an additional integer- valued state variable Q representing 
the number of exercise rights used at any considered time. 

However, our method works and the numerical results that we obtain are consistent 
with the theoretical convergence rate given in Theorem [1] One can expect that our 
method would work better on less degenerate problems. 

Remark 6 (Dealing with more exercise rights) . The computational time of our method 
intuitively increases linearly with the number of exercise rights nmax- Indeed, at each 
time step of the backward induction procedure, the number of conditional expectation 
estimations is proportional to rimax- Besides, when multiplying by 2 the number of 
exercise rights, it would require, at least, a double number of Monte Carlo samples for 
a same accuracy of the computation of conditional expectation estimators. 

Let us mention that the computational time of the benchmark method using it- 
eration increases linearly as function of the the maximal number of exercise rights as 
well. 
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